clear all; clc;close all;
%Updated: 05/07/2018
cd ..
dataA0 = load('output//mcdistA.txt');  
dataB0 = load('output_baselineX//mcdistA.txt');  

maxdata=10000;

dataA = dataA0(1:maxdata);  
dataB = dataB0(1:maxdata);  

a = load('output//agrid2.txt');  
b = load('output//bgrid2.txt');  

na=rows(a);
nb=rows(b);

[aa,bb]=meshgrid(a,b);
    
for i_a=1:na;
for i_b=1:nb;
    distAaccess(i_b,i_a) = sum(dataA== i_a +(i_b-1)*na)/maxdata;
    distBaccess(i_b,i_a) = sum(dataB== i_a +(i_b-1)*na)/maxdata;
end;
end;

figure;surf(aa,bb,distAaccess);title('model A');
figure;surf(aa,bb,distBaccess);title('model B');

sum1=sum(sum(distAaccess,2),1) 
sum2=sum(sum(distBaccess,2),1)

for i_a=1:na;
    distAaut(i_a) = sum(dataA== na*nb+i_a)/maxdata;
    distBaut(i_a) = sum(dataB== na*nb+i_a)/maxdata;
end;

figure;plot(a,distAaut);title('model A');
figure;plot(a,distBaut);title('model B');

sum1aut=sum(distAaut) 
sum2aut=sum(distBaut)
